SUBROUTINE setup_hf
  USE constants
  USE single_particle_orbits
  USE wave_functions
  USE configurations
  IMPLICIT NONE
  INTEGER :: number_orbits
  REAL(DP), ALLOCATABLE, DIMENSION(:,:)  :: coeffs  

  ALLOCATE(coeffs(all_orbit%total_orbits,all_orbit%total_orbits))

  coeffs = 0.0_dp
 
  CALL brueckner_hartree_fock(coeffs,all_orbit%total_orbits)


END


SUBROUTINE brueckner_hartree_fock(coeffs, ncoeffs)
  USE single_particle_orbits
  USE configurations
  USE constants
  USE wave_functions
  IMPLICIT NONE
  TYPE (configuration_descriptor) :: sp_configs
  INTEGER, INTENT(IN) :: ncoeffs
  REAL(DP), DIMENSION(ncoeffs,ncoeffs), INTENT(INOUT)  :: coeffs  
  INTEGER :: a, c, ket , bra, nrot, hole_states, hf_iter, max_hfiter
  INTEGER :: la, lc, na, nc, k, i,j
  INTEGER :: n_min
  REAL(DP), ALLOCATABLE, DIMENSION(:,:) :: hf_vect
  REAL(DP), ALLOCATABLE, DIMENSION(:,:) :: hf_mtx
  REAL(DP), ALLOCATABLE, DIMENSION(:,:) :: hf_kin
 
  REAL(DP), ALLOCATABLE, DIMENSION(:) :: hf_eigen
  REAL(DP), DIMENSION(n_rel) ::  sum_wf
  REAL(DP) :: e_kin, hf, sigma, total_e
  REAL(DP), ALLOCATABLE, DIMENSION(:) :: int_factor, kin_energy
  ALLOCATE (int_factor(n_rel), kin_energy(n_rel) )
  int_factor(:)=wra(:)*ra(:)*ra(:)

  !     initialize bhf harmonic oscillator sp wave function
  bhf_hol=0.0_dp; max_hfiter = 100
  ALLOCATE (hf_mtx(all_orbit%total_orbits,all_orbit%total_orbits))
  ALLOCATE (hf_vect(all_orbit%total_orbits,all_orbit%total_orbits))
  ALLOCATE (hf_kin(all_orbit%total_orbits,all_orbit%total_orbits))
  ALLOCATE (hf_eigen(all_orbit%total_orbits))
  hf_vect = 0.0_dp
  !  for first iteration coeffs has only diagonal values equal 1
  DO bra = 1, all_orbit%total_orbits
     coeffs(bra,bra) = 1.0_dp
  ENDDO
  hf_iter = 0; sigma = 1.0_dp



      hf_kin=0.0_dp
     DO bra = 1, all_orbit%total_orbits
        a = bra
        DO ket= bra, all_orbit%total_orbits
           c = ket
           IF(all_orbit%ll(a) /= all_orbit%ll(c)) CYCLE
           IF ( all_orbit%jj(a) /= all_orbit%jj(c)) CYCLE
           IF(all_orbit%itzp(a) /= all_orbit%itzp(c) ) CYCLE
           hf = 0.0_dp               


           e_kin = 0.0_dp

           n_min=min(all_orbit%nn(a),all_orbit%nn(c))+1
           if(a == c ) then
                 e_kin= (2.0_dp*n_min+all_orbit%ll(a)-0.5)*hbar_omega/2.0_dp
           elseif (abs(all_orbit%nn(a)-all_orbit%nn(c)) == 1) then
                 e_kin= sqrt(n_min*(n_min+all_orbit%ll(a)+0.5_dp))*hbar_omega/2.0_dp
           else
                 e_kin=0.0_dp
           endif
          e_kin= e_kin*real(mass_nucleus-1)/mass_nucleus
      hf_kin(a,c)=e_kin;
        ENDDO
     ENDDO


  
  DO WHILE ((hf_iter <= max_hfiter).AND.(ABS(sigma) > 1E-8))
     !     set up bhf matrix to be diagonalized
     hf_mtx = 0.0_dp; hf_vect = 0.0_dp
     DO bra = 1, all_orbit%total_orbits
        a = bra
        DO ket= bra, all_orbit%total_orbits
           c = ket
           IF(all_orbit%ll(a) /= all_orbit%ll(c)) CYCLE
           IF ( all_orbit%jj(a) /= all_orbit%jj(c)) CYCLE
           IF(all_orbit%itzp(a) /= all_orbit%itzp(c) ) CYCLE
           hf = 0.0_dp               


           CALL  diagram_HF(a, c, coeffs, hf_iter, hf, ncoeffs)
           !  compute the kinetic energy or unperturbed one-body H
           e_kin = 0.0_dp

           n_min=min(all_orbit%nn(a),all_orbit%nn(c))+1
           if(a == c ) then
                 e_kin= (2.0_dp*n_min+all_orbit%ll(a)-0.5)*hbar_omega/2.0_dp
           elseif (abs(all_orbit%nn(a)-all_orbit%nn(c)) == 1) then
                 e_kin= sqrt(n_min*(n_min+all_orbit%ll(a)+0.5_dp))*hbar_omega/2.0_dp
           else
                 e_kin=0.0_dp
           endif
          e_kin= e_kin*real(mass_nucleus-1)/mass_nucleus
           hf_mtx(bra,ket) = hf+ e_kin
           hf_mtx(ket,bra) = hf_mtx(bra,ket)
        ENDDO
     ENDDO
     !     obtain the BHF coefficients 
     CALL matrix_diag(hf_mtx,all_orbit%total_orbits, all_orbit%total_orbits,&
          hf_eigen, hf_vect,nrot)
          
          
     !     set up the new bhf harmonic  oscillator wave function
     !     Note memory stride
     sigma = 0.0_dp
     DO bra = 1, all_orbit%total_orbits
        a = bra
        la = all_orbit%ll(a)
        na = all_orbit%nn(a)
        sum_wf = 0.0_dp
        DO ket= 1, all_orbit%total_orbits
           c = ket
           lc = all_orbit%ll(c)
           nc = all_orbit%nn(c)
           sum_wf(:)=sum_wf(:)+hf_vect(ket,bra)*hol(:,lc,nc)
           coeffs(a,c) = hf_vect(bra,ket)
        ENDDO
        bhf_hol(:,a)=sum_wf(:)
        ! set up of new single-particle energy 
        sigma = sigma +ABS(all_orbit%e(a) - hf_eigen(bra))
        all_orbit%e(a) = hf_eigen(bra)                    
     ENDDO

     total_e=0.0_dp
     DO i=1,all_orbit%total_orbits
        IF(all_orbit%orbit_status(i) /= 'hole') CYCLE
        total_e= total_e + (all_orbit%e(i)+hf_kin(i,i))*(all_orbit%jj(i)+1)
    END DO
     total_e =total_e /2.0_dp
     
     sigma = sigma/all_orbit%total_orbits
     WRITE(6,*) 'iteration nr and sigma, e_total', hf_iter, sigma, total_e
     DO i=1, all_orbit%total_orbits
        IF (all_orbit%orbit_status(i) /= 'hole') CYCLE         
        WRITE(6,'(7HNumber:,6(I3,2X),2X,E20.10)') i, all_orbit%nn(i), all_orbit%ll(i), &
             all_orbit%jj(i), &
             all_orbit%nshell(i), all_orbit%itzp(i), all_orbit%e(i)
     ENDDO
     hf_iter = hf_iter+1
  ENDDO
  
  
  
  DO i=1, all_orbit%total_orbits
     IF ( keep_originalenergies == 'no') THEN 
        IF (all_orbit%model_space(i) == 'inside')  THEN
           all_orbit%evalence(i) = all_orbit%e(i)
        ELSE
           all_orbit%evalence(i) = 0.0_dp
        ENDIF
     ELSEIF ( keep_originalenergies == 'yes') THEN
        all_orbit%e(i) = all_orbit%e_original(i)
     ENDIF
     WRITE(8,'(7HNumber:,5(I4,2X),2X,E13.6,2X,E13.6,2X,A10,2X,A10)') i, all_orbit%nn(i), all_orbit%ll(i), &
          all_orbit%jj(i), all_orbit%itzp(i), all_orbit%e(i), all_orbit%evalence(i), &
          all_orbit%orbit_status(i), all_orbit%model_space(i)
  ENDDO
  DEALLOCATE ( hf_mtx)
  DEALLOCATE ( hf_eigen)
  DEALLOCATE ( hf_vect)
  !     new harmonic oscillator wave function
  wave_function=bhf_hol
  DEALLOCATE (int_factor, kin_energy )

END SUBROUTINE brueckner_hartree_fock



SUBROUTINE diagram_HF(a,c,coeffs, iteration,onebody_diagram_HF, ncoeffs)
  USE constants
  USE single_particle_orbits
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: a, c, iteration, ncoeffs
  INTEGER :: j_min, j_max, jph, h, h1, h2, i
  REAL(DP), DIMENSION(ncoeffs,ncoeffs), INTENT(IN)  :: coeffs  
  REAL(DP) :: val, ang_mom_factor, w
  REAL(DP), INTENT(INOUT) :: onebody_diagram_HF
  REAL(DP), DIMENSION(n_startenergy_g) :: ans

  onebody_diagram_HF=0.0_dp
  DO h=1, all_orbit%total_orbits
     IF (all_orbit%orbit_status(h) /= 'hole') CYCLE         
     j_min=ABS((all_orbit%jj(a)-all_orbit%jj(h))/2)
     j_max=(all_orbit%jj(a)+all_orbit%jj(h))/2
     DO h1=1, all_orbit%total_orbits
        IF  (all_orbit%jj(h) /= all_orbit%jj(h1)) CYCLE
        IF  (all_orbit%ll(h) /= all_orbit%ll(h1)) CYCLE
        IF  (all_orbit%itzp(h) /= all_orbit%itzp(h1)) CYCLE
        DO h2=1, all_orbit%total_orbits
           IF  (all_orbit%jj(h) /= all_orbit%jj(h2)) CYCLE
           IF  (all_orbit%ll(h) /= all_orbit%ll(h2)) CYCLE
           IF  (all_orbit%itzp(h) /= all_orbit%itzp(h2)) CYCLE
              DO jph=j_min,j_max                                                                    
                 ang_mom_factor=(2.*jph+1.)/(all_orbit%jj(a)+1.)                                    
                 CALL pphhmtx(a,h1,c,h2,jph,ans); IF ( ans(1) == 0.0_dp ) CYCLE                     
                 onebody_diagram_HF=onebody_diagram_HF+ans(1)*ang_mom_factor*coeffs(h1,h)*coeffs(h2,h)
              ENDDO
        ENDDO
     ENDDO
  ENDDO

END  SUBROUTINE diagram_HF



!
!     This function returns the matrix for V
!
SUBROUTINE pphhmtx(ja,jb,jc,jd,jt,gmtpn)
  USE stored_bare_interaction
  USE single_particle_orbits
  USE constants
  IMPLICIT NONE
  LOGICAL TRIAG
  INTEGER, INTENT(IN) :: ja,jb,jc,jd,jt
  INTEGER :: iph
  REAL(DP), INTENT(INOUT), DIMENSION(n_startenergy_g) :: gmtpn
  REAL(DP) :: dij, delta, xxx

  gmtpn=0.
  IF(2*all_orbit%nn(ja)+all_orbit%ll(ja) +2*all_orbit%nn(jb)+all_orbit%ll(jb) > nlmax ) RETURN
  IF(2*all_orbit%nn(jc)+all_orbit%ll(jc) +2*all_orbit%nn(jd)+all_orbit%ll(jd) > nlmax ) RETURN
  IF((all_orbit%itzp(ja)+all_orbit%itzp(jb)) /= &
       (all_orbit%itzp(jc)+all_orbit%itzp(jd))) RETURN
  IF((-1)**(all_orbit%ll(ja)+all_orbit%ll(jb)) /=  &
       (-1)**(all_orbit%ll(jc)+all_orbit%ll(jd))) RETURN
  IF((ja == jb).AND.(MOD(jt,2)/=0)) RETURN       
  IF((jc == jd).AND.(MOD(jt,2)/=0)) RETURN       
  IF(triag(all_orbit%jj(ja),all_orbit%jj(jb),2*jt)) RETURN
  IF(triag(all_orbit%jj(jc),all_orbit%jj(jd),2*jt)) RETURN
  CALL mtx_elements(ja,jb,jc,jd,jt,gmtpn)
  !  Test purpose, leave as is
  !  IF( (all_orbit%itzp(ja)+all_orbit%itzp(jb) == 0 ) ) THEN
  !     gmtpn = delta(ja,jc)*delta(jb,jd)
  !  ELSE
  !     gmtpn = (delta(ja,jc)*delta(jb,jd)-delta(ja,jd)*delta(jb,jc)*iph((all_orbit%jj(ja)+all_orbit%jj(jb))/2-jt))/dij(ja,jb)/dij(jc,jd)
  !  ENDIF

  !  gmtpn = gmtpn/dij(ja,jb)/dij(jc,jd)

END SUBROUTINE pphhmtx
!
!     Calculates the crosscoupled matrix element type 1
!
SUBROUTINE cross_coupled_mtxel1(ja,jb,jc,jd,jtot,cross)
  USE single_particle_orbits
  USE constants
  USE ang_mom_functions
  IMPLICIT NONE 
  REAL(DP),  DIMENSION (n_startenergy_g) :: ans
  INTEGER :: ja,jb,jc,jd,jtot,jbra_min, jket_min, jbra_max ,jket_max, &
       jt, j_min, j_max, iph
  REAL(DP) :: fnorm, ang_mom_factor
  REAL(DP), DIMENSION(n_startenergy_g), INTENT(OUT) :: cross

  cross=0.
  fnorm=SQRT(2.*jtot+1.)*iph((all_orbit%jj(ja)+all_orbit%jj(jd))/2+jtot)
  IF(iph(all_orbit%ll(ja)+all_orbit%ll(jb)) & 
       /=iph(all_orbit%ll(jc)+all_orbit%ll(jd))) RETURN
  IF((all_orbit%itzp(ja)+all_orbit%itzp(jb)) /= &
       (all_orbit%itzp(jc)+all_orbit%itzp(jd))) RETURN
  jbra_min=ABS(all_orbit%jj(ja)-all_orbit%jj(jb))
  jbra_max=all_orbit%jj(ja)+all_orbit%jj(jb)
  jket_min=ABS(all_orbit%jj(jc)-all_orbit%jj(jd))
  jket_max=all_orbit%jj(jc)+all_orbit%jj(jd)
  j_max=MIN(jbra_max,jket_max)/2
  j_min=MAX(jbra_min,jket_min)/2
  IF(j_min > j_max) RETURN
  DO jt=j_min,j_max
     ang_mom_factor=sjs(all_orbit%jj(jc),all_orbit%jj(ja),2*jtot, &
          all_orbit%jj(jb),all_orbit%jj(jd),2*jt)* &
          (2.*jt+1.)*fnorm*iph(jt)
     CALL pphhmtx(ja,jb,jc,jd,jt,ans)
     cross=cross+ans*ang_mom_factor
  ENDDO

END SUBROUTINE cross_coupled_mtxel1
!
!     Calculates the crosscoupled matrix element type 2
!
SUBROUTINE cross_coupled_mtxel2(ja,jb,jc,jd,jtot,cross)
  USE single_particle_orbits
  USE constants
  USE ang_mom_functions
  IMPLICIT NONE 
  REAL(DP), DIMENSION (n_startenergy_g) :: ans
  INTEGER :: ja,jb,jc,jd,jtot, jbra_min, jket_min, jbra_max ,jket_max, &
       jt, j_min, j_max, iph
  REAL(DP) :: fnorm, ang_mom_factor
  REAL(DP), DIMENSION(n_startenergy_g), INTENT(OUT) :: cross

  cross=0.
  fnorm=SQRT(2.*jtot+1.)*iph((all_orbit%jj(ja)+all_orbit%jj(jd))/2+ &
       jtot+all_orbit%jj(jb))
  IF(iph(all_orbit%ll(ja)+all_orbit%ll(jb)) & 
       /=iph(all_orbit%ll(jc)+all_orbit%ll(jd))) RETURN
  IF((all_orbit%itzp(ja)+all_orbit%itzp(jb)) /= &
       (all_orbit%itzp(jc)+all_orbit%itzp(jd))) RETURN
  jbra_min=ABS(all_orbit%jj(ja)-all_orbit%jj(jb))
  jbra_max=all_orbit%jj(ja)+all_orbit%jj(jb)
  jket_min=ABS(all_orbit%jj(jc)-all_orbit%jj(jd))
  jket_max=all_orbit%jj(jc)+all_orbit%jj(jd)
  j_max=MIN(jbra_max,jket_max)/2
  j_min=MAX(jbra_min,jket_min)/2
  IF(j_min > j_max) RETURN
  DO jt=j_min,j_max
     ang_mom_factor=sjs(all_orbit%jj(jc),all_orbit%jj(jb),2*jtot, &
          all_orbit%jj(ja),all_orbit%jj(jd),2*jt)* &
          (2.*jt+1.)*fnorm
     CALL pphhmtx(ja,jb,jc,jd,jt,ans)
     cross=cross+ans*ang_mom_factor
  ENDDO

END SUBROUTINE cross_coupled_mtxel2
!
!Fast diagonalization of real(kind=8) symmetric matrices
!
SUBROUTINE diag_exact(cvec, ceig, n, h)
  USE constants
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: n
  COMPLEX*16, DIMENSION(n,n), INTENT(IN) :: h
  INTEGER :: k, np, p, i,j, i1, kvec, lwork, option, il ,iu, info
  REAL(DP) :: a(n*(n+1)), work(30*n), thresh, vl, vu
  LOGICAL :: flag(n)  
  COMPLEX*16 :: cvu, cvl
  COMPLEX*16, INTENT(INOUT) :: ceig(n), cvec(n,n)

  cvec = 0.0_dp; ceig = 0.0_dp
  lwork = 30*n
  thresh = 30.0
  kvec = n
  option = 4
  info = 1
  i1 = 0
  DO i =  1, n
     DO j =  1, i
        i1 = i1 + 1
        a(i1) = DBLE(h(j,i))
        a(i1+n*(n+1)/2) = AIMAG(h(j,i))
     ENDDO
  ENDDO
  CALL cs(n,a,ceig,kvec,cvec,lwork,work,thresh, option,il,iu,cvl,cvu,flag,info)

END SUBROUTINE diag_exact
!
!Fast diagonalization of real(kind=8) symmetric matrices
!
SUBROUTINE diag_exactzgeev(cvec, ceig, n, h)
  USE constants
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: n
  COMPLEX*16, DIMENSION(n,n), INTENT(IN) :: h
  COMPLEX*16, DIMENSION(n,n), INTENT(INOUT) :: cvec
  COMPLEX*16, DIMENSION(n), INTENT(INOUT) :: ceig
  REAL(DP), DIMENSION(2*n) :: rwork
  COMPLEX*16, DIMENSION(n,n) :: vl
  COMPLEX*16, DIMENSION(10000) :: work1
  INTEGER :: i, lda, ldb, ldvl, ldvr, lwork, info
  CHARACTER*1 :: jobvl, jobvr, balanc, sense

  balanc = 'n';  jobvl = 'n' ;  jobvr = 'v';  sense = 'n';  lda = n
  ldvl = 1;  ldvr = n;  lwork = 10000
  cvec = 0.0_dp; ceig = 0.0_dp
  CALL zgeev( jobvl, jobvr, n, h, lda, ceig, vl, ldvl, cvec, ldvr, &
       work1, lwork, rwork, info )

END SUBROUTINE diag_exactzgeev
!
! Get eigenvalues for large space, eigenvectors for model space and excluded space
!
SUBROUTINE eigenvalues_large(cvec_pp,cvec_qp,ceig_p,hamilton,n,np)
  USE single_particle_orbits
  USE configurations
  USE constants
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: n, np
  INTEGER :: k1, i, nqq
  COMPLEX*16, ALLOCATABLE, DIMENSION(:,:) :: cvec
  COMPLEX*16, ALLOCATABLE, DIMENSION(:) :: ceig
  COMPLEX*16, INTENT(INOUT)  :: ceig_p(np)
  COMPLEX*16, INTENT(INOUT) :: cvec_pp(np,np)
  COMPLEX*16, INTENT(INOUT) :: cvec_qp(n-np, np) 
  COMPLEX*16, INTENT(IN) :: hamilton(n,n)
  REAL(DP) :: cvec_max
  REAL(DP), ALLOCATABLE, DIMENSION(:) :: temp

  ! diagonalize 2p-effective shell model hamiltonian, this takes a symmetric 
  ! matrix as input.
  ALLOCATE(cvec(n,n)); ALLOCATE(ceig(n))
  cvec = 0.0_dp; ceig = 0.0_dp
  CALL diag_exact(cvec, ceig, n, hamilton)
  cvec_pp = 0.0_dp; cvec_qp = 0.0_dp
  ! Set up eigenvalues and eigenvectors for model space and excluded space
  DO k1=1, np
     ! loop over all model space coefficients of exact eigenvectors |k>
     nqq = 0
     DO i = 1, n
        IF ( i <= np ) THEN
           cvec_pp(i,k1) = cvec(i,k1) 
           ceig_p(k1) = ceig(k1) 
        ELSEIF ( i > np ) THEN
           nqq = nqq + 1
           cvec_qp(i-np,k1) = cvec(i,k1) 
        ENDIF
     ENDDO
  ENDDO
  DEALLOCATE(ceig); DEALLOCATE(cvec)

END SUBROUTINE eigenvalues_large
